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FOREWORD 

This report presents the analysis and the computer programs 
developed for computation of viscous shock layer flowfield 
surrounding the nose of a shuttle orbiter during its reentry 
Part I describes the problem formulation and the numerical 
procedures used to solve the basic set of equations, and the 
results of flowfield properties at several trajectory points 
ranging from the high altitude rarefied region to the low 
altitude boundary-layer region. Part II of this report 
describes the structure of the computer programs and the 
experiences gained in utilizing these programs. A user's 
input guide is also included along with a complete listing 
of programs . 
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PART I 

ANALYSIS AND RESULTS FOR A SHUTTLE ORBITER 
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SUMMARY 

This part of the report presents a numerical analysis of 
the viscous flow field, with and without finite-rate chemical 
reactions, in the nose region of a shuttle orbiter under a 
wide range of free stream conditions. One purpose of this 
study was to develop a unified calculation procedure that 
will provide a starting solution having detailed profiles 
of flow properties for subsequent flow field computations 
beyond the nose. Therefore, the generalized curvilinear 
coordinate system was used and the fluid-dynamic equations 
were cast in conservative form. Thus, several special 
coordinate systems can be chosen in the computation and 
the shock can be treated as either a sharp discontinuity 
or a thick layer. The second objective of this study was 
to investigate the flow field characteristics that are 
encountered during the orbiter descent. The effect of 
transport properties of the air mixture, the surface 
catalyticity , and the wall temperature on the flow field 
was studied extensively at several trajectory points for 
which the chemical nonequilibrium phenomena are predominant. 
The last objective was to analyze the flowfield in terms 
of the heat transfer and friction coefficients and to 
compare the results with available solutions. Representive 
trajectory points were selected for calculations using 
the frozen, finite-rate, and equilibrium gas models. The 
numerical solutions obtained are considered to be sufficiently 
accurate for the aforementioned objectives due to the use 
of exact equations, and the coordinate transformation which 
provides a better resolution of flow properties in the 
vicinity of a wall. Attempts were also made to improve the 
efficiency of the time-marching finite-difference technique 
which was used to solve the flow equation in the present 
analysis . 
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L 
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NOMENCLATURE 

frozen speed of sound 

temperature coefficients for 
polynomial equations of thermodynamic 
functions, Eq (16) 

body configuration 

mass fraction 

friction coefficient, Eq (24) 

heat transfer coefficient (Stanton 
number) , Eq (25) 

heat capacity at constant pressure 

binary diffusion coefficient 

multicomponent diffusion coefficient 

specific internal energy, Eq (6) 

Gibbs free energy Eq (16) 
specific enthalpy 

metric coefficients for the orthogonal 
coordinate system 

molar enthalpy, Eq (16) 

mass diffusion flux, Eq (9) 

thermal conductivity, Eq (18) and (19) 

rate constants of forward reaction, 
backward reaction 

Boltzmann constant 

equilibrium constant for mass concentrations 

total number of species 

Lewis number for a multicomponent 
mixture 

molecular weight or Mach number 
molecular weight of the mixture 
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00 
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NOMENCLATURE 

normal to wall or shock 
pressure, Eq (5) 

Prandtl number 
heat flux, Eq (8) 

internal partition function 
universal gas constant 
Reynolds number 
nose radius 

tangent to wall or shock 

specific entropy or shock configuration 

time coordinate 

temperature 

velocity components 

free stream speed 

mole fraction 

space coordinates 

flow quantity vector 

parameter in stretching coordinate 


e 

6 


Y 

0) 

n 

p 


O U/m) 

u ik 


e + % (u 2 +v 2 ) 

distance between shock and body along 
£-axis 

ratio of specific heats 
chemical rate of production 
stress tensor, Eq (7) 
density 

defined in Eq (7) 
viscosity, Eq (17) 

collision cross section 
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K, n 

Subscripts 

i/ j 

i ,m 

00 

T 

w 

Superscripts 


NOMENCLATURE 

N 

space coordinates in orthogonal 
coordinate system 


space dimensions 
£-th, m-th species 
free stream 
total condition 
wall 


quantities normalized by free stream 
conditions 

flow variables in the transformed 
computational plane 
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1.0 INTRODUCTION 

This document is the final report summarizing the work 
accomplished for Tasks G and K of project 3782. The primary 
purpose of these tasks was to study and to develop a program 
for analyzing the chemical nonequilibrium, viscous flowfield 
in the nose region of a shuttle orbiter. As discussed in 
Ref. 1 the full set of Navier-Stokes (N-S) equations is 
used and the chemical kinetic equations are coupled to the 
N-S equations to achieve a higher accuracy in the flowfield 
analysis. 

1.1 Objectives of the Study 

There are three particular areas which were pursued in the 
course of this study in order to establish an efficient 
and accurate calculation scheme. The first area of interest 
was to formulate the governing equations in an unsteady 
conservative form using generalized curvilinear coordinates . 
The conservative form of equations possesses not only mathe- 
matical simplicity, but the capability to determine the 
imbedded shocks in the flow. This feature is needed to 
compute the flowfield at high altitudes where the shock is 
no longer a thin layer. The equations being written in 
curvilinear coordinates also faciliates the flowfield com- 
putation for several particular coordinate systems that can 
be chosen to define a starting line for subsequent super- 
sonic flow calculation downstream of the nose. 

The second area of concern was to develop a self-contained 
procedure for the calculations of thermodynamic and transport 
properties of the air mixture. Existing procedures are 
either limited to a certain range of temperatures and/or 
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involve approximations which have not been justified in 
the calculation of transport properties. This study there- 
fore uses a more satisfactory procedure that is at least 
as accurate as the flow model itself. 

The last area of interest was to calculate the blunt body 
flow at several representative flight conditions to cover 
the entire flight trajectory. A wide spectrum of flow 
characteristics exist in the descent trajectory; namely, 
the classical boundary layer regime at low altitude, the 
nonequilibrium shock and boundary layer interactive region 
in the middle portion of the trajectory, and the rarefied 
flow environment at high altitude. Since the present 
analysis is intended to provide the complete flow solution 

around the orbiter nose, therefore, the flow field solutions are 
obtained within the scope of a continuum flow model. 
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1.2 Numerical Analyses for Viscous Reacting Flow fields 

Because of past interest in the reentry technology and the 
associated experimental study of reacting gas flows, numerous 
analyses are available for analyzing the reacting viscous 
flow around a blunt body. Solutions have been obtained for 
both high and low Reynolds number limits because the de- 
parture from chemical equilibrium can significantly effect 
the flow observables and the skin heat transfer coefficient 
over a broad range of altitudes and speeds. The three most 
popular numerical methods developed in the last decade are 
the finite-difference method , ^ 2 9 ^ difference-differential 
method, ( 10 J , and the method of integral relations (11 13) . 
These methods have been used quite extensively in the in- 
vestigation of non-reacting viscous problem including the 
boundary-layer and the thin shock layer, and considerable 
successes have been obtained. However, in dealing with 
the reacting flow problem where the chemical nonequilibrium 
processes couple directly with the fluid-dynamic equations 
these methods are not as successful. The basic difficulty 
lies in the fact that certain assumptions of the flow must 
be met, or some input data must be given in order to carry 
out the f lowfield analysis. For boundary-layer analyses 
the edge condtions for all dependent variables should be 
specified before one can use any of the three methods. 

The boundary layer edge location can not be simply defined 
by the inviscid nonequilibrium calculation as it can for 
the non-reacting case, since the swallowing of inviscid flow 
has to be considered ' . To obtain an accurate result 

from the nonequilibrium boundary layer analysis, several 
iterations are required between the outer inviscid flow 
calculation and the boundary layer calculation. Thus, the 


10 



TR2007 


modeling of flow field by an outer inviscid layer and an 
inner boundary layer, using the three highly developed 
numerical methods becomes less attractive. 

The current analyses of the thin viscous shock layers are 
not as sophisticated as that of the boundary layer. Many 
analyses only provide solutions along the stagnation stream 
line, (4,7,9) , an( j some analyses which yield solutions for 
the shock layer were developed using many simplifying 
assumptions for defining the transport properties, and gas 
models (S" 9 ' 12 ' 13 ). i n addition to the basic assumption 
that the thickness of the shock layer is very thin compared 
to the nose radius, the shock wave itself has to be treated 
as a mathematical discontinuity. At high altitude this 
assumption gradually becomes less justified. Velocity slip 
is then introduced to the usual Rankine-Hugoniot relations 
with the hope that the downstream flow properties can be 
determined with acceptable accuracy. This scheme, known 
as the "two- layer" model in the literature (6) , has provided 
valuable results up to the point where the shock wave and 
the boundary layer merges with each other and as long as there 
is no departure of chemical equilibrium in the flow. How- 
ever, it appears that the chemical nonequilibrium processes 
would incur a certain amount of ambiguity in the shock wave 
calculations, as has been demonstrated in the various 
nonequilibrium shock layer analysis published. 

On the basis of this proceeding discussion, we have come 
to the understanding that the finite-difference method 
and the method of integral relations have been developed 
mainly for the boundary- layer type equations. The necessary 
conditions for such analyses to be valid are that the thick- 
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ness of the boundary layer or shock layer must be thin 
compared to the nose radius and that boundary .conditions be 
specified on the edge of the layer. These requirements 
reduce the applications of these methods to a rather small 
portion on the entire orbiter trajectory. Therefore, a 
satisfactory theory that can analyze the flow field for the 
flight altitudes between 200 and 300K FT, where the orbiter 
will experience significant aerodynamic heating, is still 
in demand . 

The time-marching finite-difference method used in the 
present study has received considerable attention in recent 
years. It appeals to the flow field analyst mainly because 
the exact governing equations can be used, and the accuracy 
of the solution is dependent on the mesh size only. The 
essential concept of this method is to simulate the flow 
field development from a given set of initial conditions 
until the flow settles down to its steady state. Although 
the steady solution is what one seeks , the introduction of 
the unsteady term in the equations is necessary from the 
mathematical point of view, because it changes the parabolic 
or the elliptical type of equation into a hyperbolic type for 
which a powerful numerical method is available. This method 
has been used by many investigators to solve inviscid flow 
problem of practical interest, including cases of high angle 
of attack, finite-rate chemical reaction and radiative heat 
transfer. But there exist few practical applications for 
the viscous flow problem. Most of the applications have 
been two-dimensional problems with a simple gas model, 
presumably due to the requirement of long computer time 
and the fact that the previous viscous, nonequilibrium flow 
analyses are quite adequate for designing many re-entry vehicles. 
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The presently, more strigent, design of a reusable orbiter 
requires more exact flow field analyses, that can only be 
achieved by means of the time-marching finite-difference 
method. 

A partial list of recent publications on the subject of 
viscous flow problems and their approximate range of 
validity is given in Table 1. It is not intended to com- 
prise a complete bibliography of significant research in 
this area, but to indicate the state-of-art of the viscous 
flow analyses under these groups; namely, the reacting 
multicomponent and binary mixtures, and the non-reacting 
gas. It is seen that although the previous analyses can 
be used to cover a wide range of flight conditions, such a 
complete solution can only be obtained for the stagnation 
region. 

A comparison of the afore-discussed numerical method's can be 
made on the basis of their applications to the multicomponent 
nonequilibrium blunt body flow problem. The flow characteris- 
tics of interest are the shock location, and the properties 
inside the shock layer and along the body downstream of the 
stagnation point. The calculation of the stagnation flow 
is not included in this comparison. It is seen in Table 2 
that the analysis based on the time-marching method may 
represent a unified numerical approach and provide more 
satisfactory results because of the coupling of chemical 
kinetic equations and the Navier-Stokes equations, and the 
use of the rigorous theory on the transport properties. The 
ease of application and the cost of computer time are also 
listed in approximate terms, to aid one in making an 
evaluation of the three methods. 
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2.0 GOVERNING EQUATIONS AND BOUNDARY CONDITIONS 

2 . 1 Assumptions and Governing Equations 

The present flow field analysis is based on the following 
major assumptions: 

1) Navier-Stokes model and non-slip velocity and no- 
jump temperature on the wall; 

2) negligible radiation transfer; 

3) flow in both vibrational and rotational equilibria; 

4) gas consists of a mixture of perfect gases; 

5) heat flux and mass diffusion flux are approximated 
by the Fourier and Fick laws, respectively; 

6) transport coefficients are derived for a multi- 
component mixture; 

7) negligible bulk viscosity. 

The preceeding assumptions except the first were used in 
part or wholly by previous analyses dealing with a thin 
viscous shock layer or a boundary layer. They define a 
reasonably realistic model for the flow of interest, and allow 
the solution of the model to be manageable. 
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where is the net rate of production of species 1 as a 
result of chemical reaction, e^ is the specific internal 

energy and is defined as e 1 = J c v dT. The stress tensor is 
defined as ^ 
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The heat flux vector is defined as 
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The mass diffusion vector for species 1 is defined as 


I 


1. 

x 



M 1 3 c 

= *M“^ pD lm 3x“ 
wm l 


( 9 ) 


Seven unknown scalers and vectors appear in Eqn (1-6) , 
i.e., p, u^, p, e, c^, p and T. All unknowns except T 
can be determined explicitly from the set of equations. 

T must be calculated iteratively using Eq (6) . Since a 
mixture of reacting pure gases is considered here, e^ is 
treated as a function of temperature. A polynominal 
function of T is given for e^ in Section 3.1. 

In Eqn (7-9), transport coefficients y,k and D, are 

_ lm 

determined from T, p, C^, and C3 (1 m) . The method of 

calculation of transport coefficients is discussed in 
section 3.2. 


2 . 2 Orthogonal Coordinate System 

Governing equations of Eqn (1-4) can be recast in the 
generalized orthogonal coordinates. Making use of the 
relations shown in Ref. 15 for a set of non-conservative 
form of equation, Eqn (1-4) become 


9 U 
9 1 


+ 



9G 

an 


+ H 


0 


( 10 ) 


where U, F, G and H are the column vectors that have the 
following expressions: 

U — h^h^h;} ( P f P u , p v , p £ , p c^ ) • 
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are the space coordinates. The stress and heat flux are 
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Equation (10) reduces to a simpler form along the axis of 
symmetry for an axisymmetric flow. 
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The singularity has been eliminated from Eqn (11) since h. 
no longer appears in the column vectors. 


The space coordinates and the metric coefficients for three 
^iffs^snt versions of coordinate systems in the generalized 
orthogonal coordinate system are given in Table 3 and 
illustrated in Figure 1. 


2.3 Transformation to Computational Planes 

Conventional parameters such as the Reynolds number, the 
Prandtl number, and the Lewis number can be introduced in 
Eqn (10) by the following normalization procedure. Let 
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t' = t/(Vp~ /p 00 / 1 ^) be the nondimensional independent 

variables. Substituting these quantities into Eqn (10), 
one obtains the dimensionless governing equations which have 
the same form as Eqn (10) except for the parameters appearing 
with the stress tensor, heat flux and mass diffusion vectors. 


where 
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The governing equations are then transformed to a computa- 
tional coordinate system, on which both the shock (or the outer 
surface) and the body are made boundary mesh lines of the 
computational region. Let y = 5, z = l- ~ where 6=S-B; 

the distance between the shock and body, or between the 
outer surface and the body (see figure 1) . The transformed 
equation becomes 


+ |— ( (1- z) (6.U+6 F) -G )+ 6H = 0 (12) 

at ay a z t y 

at every point except the axis of symmetry, on which the 
following form is used: 

f~ + 2^ + ((1-z) (6 t U+S y F) - G) + SH = 0 (13) 
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The stress tensor, heat flux and mass diffusion vectors in 
the computational coordinates are: 
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The conversion from the generalized orthogonal coordinate 
system to its special versions; namely, the body intrinsic, 
polar and cylindrical coordinate systems can be easily made 
with the use of Table 4. 
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The computational region is further mapped to another 
plane to allow higher space resolution near the body. 

This transformation of coordinates is desirable for an 
accurate calculation of momentum and energy transfer from the 
flow to the body. Let y and z be the coordinates of the 
second computational plane, and z = (1-exp ($z) )/ (1-exp ( 6 )) 


y = y. 
Ref. 14. 


The relations have been shown to be valuable in 


Making use of 


= -eexp(Bz) 3_ and „ 9_ 

3z 1-exp (6) 3 y 
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one can obtain the governing equations cast in the new 


computational plane as follows: 


+ i£ + H H = o 


for every where in the computational region except on the 
axis of symmetry; and 


+ 2 ^ + 1 ® + H = 0 ' (15) 

■ 9y 9z 

for the axis of symmetry. 

The variables U 3 F 3 G 3 and H and variables with overline 
are defined by: 

U = 6U, F = 6F, G = i!exp(^ } ((1 ~ Z) ^ (5 t +<S y F) ~ G) ' 

H = 6H - 6 ( (1-z) (6 t +6 y F)-G) 

U = U/ h 3 , F = F/ h 3 , G = G/hy I = 6H - 3 ( (1-z) (6 fc +6 F)-G) 

A detailed expression of Eqns (14) and (15) can be worked 
out by substituting U, F, G, H, U, F, G, H into Eqns (12) 
and (13) . 
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2 . 4 Boundary Conditions 

The velocity non-slip and temperature no-jump conditions are 
used in this analysis as the wall conditions. 

for an isothermal wall 
for an adiabatic wall 

Where e and T are related by e = c T for a highly cooled 
w w J wvw 

wall. 

The wall catalycity is only considered in two limiting forms, 
i.e. ; 


u = 0 = v 
T = T 



c. = c. (T ) for a fully catalytic wall 
X X w 

3 c • 

1=0 for a non-catalytic wall 

9n 

The boundary conditions on the outer surface of the com- 
putational region are the free stream conditions when the 
shock is to be computed inside the region. If the shock 
is used as the outer surface then the Rankine-Hugoniot 
relations are employed to determine the flow variables 
immediately behind the shock. The assumption of frozen 
composition of species across the shock is also utilized. 
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3.0 THERMODYNAMIC AND TRANSPORT PROPERTIES AND CHEMICAL 
KINETICS 

3 . 1 Method of Calculation of Thermodynamic Properties of 
Perfect Gases 

Numerous methods are presently available for the calcula- 
tion of the thermodynamic properties of species of perfect 
gases. The program of Reference 16, which was written at 
the NASA Lewis Research Center, was used as a matter of 
convenience. The difference between the method used in Ref. 

16 and others include one or more of the following: different 
forms for the partition function, different spectroscopic 
data, inclusion of excited-state data, and different heats 
of formation. The major aspects involved in the calculation 
of thermodynamic properties of the air species will be 
briefly included here, in order to present the analysis in a 
self-contained manner. More complete details can be found 
in Ref. 16. 

Equations for evaluating thermodynamic functions from the 
partition function and its first and second derivatives are 


C p_ T 2 d 2 Q ,T dQ> 2 2T dQ 5 

R Q dT 2 Q dT ; + Q dT + 2 


H T -H 0 T dQ + 5 

RT Q dT 2 


- F T -H 0 = InQ + | 


RT 


2 In M + j In T 


- 3.66511 


The internal partition function Q contains vibrational, 
rotational, and electronic contributions. The last term 
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in each of the equations is the translational contribution. 
The partition function for monatomic gases is given by 

Q - £ ex ? <-V KT) 

m=l 

where q and E are the statistical weight and electronic 
m m 

energy of the m state, respectively. The way to terminate 
the number of energy levels L is to include all energy 
levels that are less than or equal to the ionization poten- 
tial lowered by an amount KT. This cutoff method is 
temperature dependent and is used in this study. 

For diatomic gases, Q involves vibrational and rotational 
as well as electronic energy. 


Q = Q m Q m QS Q 

e V v R corr 


The quantities Q™, Q™, Q™ are the electronic, harmonic 

oscillator, and classical rotation contributions to the 
partition function, respectively. The remaining term is 
the correction term given by the modified Pennington and 
Kobe method (Ref . 29 in Ref . 16) . 


After calculating the thermodynamic data for the species, 
a least-square technique is used to fit these data into 
polynomials . The input spectroscopic data and the resultant 
coefficients of the empirical equations for thermodynamics 
functions are given in Appendix I. The following shows 
the thermodynamic functions in terms of the coefficients: 
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C /R 
P 


2 


3 m 2 


3 

+ a 4 t 

4 

+ a 5 t* 



A 4 3 

+ -j-2-T 

A 5 4 
+ —T 


(16) 

4 

5 

T 


A 3 t 2 

A 4 t 3 

As t 4 + Ag 

— A 

6“ T _ 

I2 T 

20 T 

/ 


H/RT = A 1 
F/RT = h 1 

3 . 2 Method of Calculations of Transport Property of a Mixture 

The Chapman-Enskog theory, as extended to a multi-component 

. ( 1 ?) 
mixture by Herschfelder , etc., , is used to calculate the 

transport coefficients. 

The first approximation to the coefficient of viscosity, 
for a mixture of gases, is given as 


v = 


F li 

F 12 

F 13 • 

* * F ln 

X 1 

F 12 

F 22 

F 23 * 

* ‘ F 2n 

X 2 

F 13 

F 2 3 

F 33 ' 

' * F 3n 

X 3 

• 

F ln 

F 

2n 

F 

r 3n * 

. . F 

nn 

• 

X 

n 

x. 

X. 

x _ 

. . X 

0 

1 

2 

3 

n 



(17) 


r— 1 
rH 

F 12 

F 13 • ' 

. . F. 
In 

F 12 

F 22 

F 23 * ‘ 

. . F 0 
2n 

F 13 

• 

F 

23 

• 

F 

33 * ' 

• 

. . F_ 
3n 

• 

• 

• 

F-i 

In 

• 

• 

Fn 

2n 

• 

• 

F, . . 
3n 

• 

’ * F nn 
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where the principal diagonal elements are given by 

2 
y 

F 


S ( 2 *i*k\ M A 

' n ik / (M. + M 


M, 


k^i 


i + v “ \ 3A ik 


/ 5 ~k 

T-\77— + m 7 
i , 


and the off-diagonal elements by 


'2x . x . 
_ - 1 J 


M.M. 
i D 


F ij F ji \ hj j / tt: ttttt 


n ij / (M. + M.) 

i D 


3A* . 
' ID 


- 1 ,i ^ j 


The quantity n^ k in above equations is given by 


26 


n ik* 10 = 


. / 2M.M 

- 693 y^f~ + 


V 


~¥ Trrr 

“ik 


The quantity n ii represents the viscosity coefficient for 
molecule i and may be obtained by letting k=i in the same 
equation. 

The first approximation to the coefficient of thermal con-' 
ductivity of a mixture of reacting gases contains two terms 


k = k 
An expression 


monatomic + ^internal 

f 03T k ■ 

monatomic is given as 
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'‘monatomic = 4 


11 L 12 L 13 * * 

• L ln X 1 

21 L 22 L 23 * * 

* L 2n X 2 

31 L 32 L 33 * * 

• • 

* L 3n X 3 

• • 

• • 

• • 

nl L n2 L n3 * * 

• • 

• • 

* ^nn x n 

1 x 2 x 3 * * 

• x n 0 

L 11 L 12 L 13 • 

* * L ln 

L 21 L 22 L 23 * 

* * L 2n 

L 31 L 32 L 33 * 

• • • 

• • L A 

3n 

• 

• • • 
L nl L n2 L n3 

• 

• 

L 

nn 


(18) 


where the principal diagonal elements are given by 


= V' 16x i x k(r^'i + T*i - 3M * B Ik + 4M A A ik) M A 


11 


15Rn ii 


£ 

k^i 


15R< M i. + V A ik n ik 


and the off-diagonal elements by 


L. . 
il 


L. . 
li 


16x.x.M?M? 
1 1 1 1 


isrOvl + 


Mj) 


A ^ii 


(P - 3B* . - 4A* . ) i 
\ 4 i] 11/ 


i t* j 
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An expression for k. , ..is given as 
c internal 3 


'"mixture internal 


n 6** / 5 R \ 

5 A ii ( C p. " 2 M. n ii ) 


i=l 


At . 2M . , n . . x . 

At . \ M. + M./ n • • X. 
ID l 3 'lD l 


( 19 ) 


1 + 


Finally, the first approximation to the coefficient of 
binary diffusion is given as 

3 (M. + M.) 

°ij 5 p A ij n ij 


( 20 ) 


The multicomponent diffusion coefficients are calculated by 
the following formula 


r, _ 1 „ K JX - K 11 
D ij Mj M w jlTf 


( 21 ) 


where K . . = 0 

li 


K . . = x i + 

13 *77 


M 


-1 T 


M i Jg^i D ik 


K| is the determinant of the ^ and K 1] are the minors 
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0 . 


1 = (- 1 ) 1+ ^ 


K i ■ i K, 

1 , 1-1 1 , 1+1 


K j-l,l* * K j-l,i-l K j-l,i+l 
K j+l,l* * K j+1, i-1, K j+l,i+l 


K. t 
D- l,n 

K ... 
D+l,n 


K 


n, 1 


. K 


n, i-1 


K 


n, i+1 


. . . K 


n,n 


Note that all the transport coefficients are expressed in 
terms of the quantity to facilitate computer calculations. 

Simple derivations and references based on which these 
formulae are obtained are shown in Reference (17). 


The calculation procedure for solving y and k monatomic is 

essentially the same one used in Reference 18 in which the 

equations (17) and (18) are written as a set of simultaneous 

linear algebraic equations and the Gauss-Jordan reduction 

scheme is used. In situations where the mole fraction of 

some species is zero, such as in the frozen flow, problems 

occur in solving the algebraic equations. In this case, the 

— 8 

mole fraction is set to 10 in order to avoid the round-off 
errors in the calculation and to prevent more than . negligible 
contributions to the results. 
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The binary collision cross sections and their ratios are 
needed in the calculation of transport properties. 


= n (2,2) / n (1 ' 1) 
ik / ik 


_ ( c o"^ 1 / 2 ) 4 tt( 1,3) w -r(l,l) 

B ik “ (5fi ik 4n ik )/Q ik 


There are 25 possible binary interactions for dissociating 
air having five species for which the data of cross section 

are generally available up to 10000°K. Extrapolation 
is used whenever the temperature is outside of the given range. 
In case the binary cross section is not available, a simple 
combination rule is used; i.e. 


- <A ii + \* )/2 


B ik = (B ii + \k )/2 


When the cross section is unknown for two like neutral species, 
the rigid sphere cross section is used. Since the data are 
scarce for changed species, their contribution to the transport 
properties are not accounted for, consequently the calculation 
of transport properties as described in this section is only 
satisfactory when the ionization is not significant. The 
data of cross sections used in this study are summarized in 
Appendix II. 

3 . 3 Chemical Reactions 

The calculation of the net rate of production ^ for 
species 1 can be carried out in a number of ways because 
of the number of the relations governing the chemical 
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reaction rate is more than the number of species. The 
redundancy of the chemical system is caused by two extra 
relations for the conservation of atoms or nuclei of the 
molecular species and the conservation of charge for ionized 
species in addition to the chemical reactions. The usual 
approach is to obtain the of the molecular and the 
ionized species from the chemical reactions, and the 
of the atomic species and electrons from the conservation 
relations. In dealing with a complicated chemical system 
involving a large number of species, however, it is found 
in Ref. 19 that the calculation of simply from the 
chemical reactions alone is more expedient. Since the 
conservation of atoms and charge is satisfied with a chosen 
set of chemical reactions, the unique solutions of is 
obtained. 


In the present analysis six different types of reactions 
can be considered. They are listed as follows: 


I 

A+B^ C+D 

II 

A+B+(M)i C+(M) 

III 

A+B^ C+D+E 

IV 

A+B^C 

V 

A+(M) B+C t (M) 

VI 

A+B+C ^D+E 

A, B , C, D and 

E refer to the reacting species and M denotes 


the third body. The production and dimunition rate of 
species involved in the reactions are respectively (w 1 ) f 
and . 
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("l’b 


I. 

k fP 2 [A] [B] 

k b p 2 [C] [D] 


II. 

k f p 3 (A3 (B]/M w 

k b p 2 [C]/M w /RT 

III. 

k f p 2 [A] [B] 

k. p 3 [C] [D] [E] RT 
b 

• 

> 

M 

k f p 2 [A] [B] 

k b p [C] /RT 


V. 

k f p 2 [A]/M w 

k b p 3 [B] [C]RT/M w 

VI. 

k b p 3 [A] [B] [C] RT 

k f p 2 [D] [E] 


where 

k = k./K , K is 
b f c c 

the equilibrium constant 


k f = A T n exp(- e/ T) 


r . 3 ~\ 

cm 



f 

mole sec 



the coefficients A, n , and e are given Appendix III, 

where k f is usually given in literature in the cgs-Kelvin 

unit. In the formulas for calculations of (cu.)^ and (u).), 

it Id 

[ ] indicates the mass concentration per mole of the species 

3 

and R = 82.07835 atm-cm /°K/mole. The net rate of production 
is obtained by the following relation to reduce the round-off 
error 

N ' N v 

W 1 = (X1 U l ) f,r ■ Lj U l ) b,r 
“ r=1 

where r refers to the number of chemical 


) M x (22) 

P 

reactions . 
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4.0 TIME-MARCHING FINITE -DIFFERENCE METHOD 

4 . 1 Predictor-Corrector Technique 

Among several versions of the Lax-Wendroff second-order 

difference scheme for solving an initial-value problem, the 

20 

predictor-corrector technique proposed by MacCormack is 
probably the most widely used. This technique not only 
yields better accuracy, but is easier to use than other 
versions, especially for mesh points located on the compu- 
tational boundaries. The essence of this technique can be 
described as follows: 


k+1 

x = 

n,m 

x k + 

n,m 

, 3x, k 
9 1 n ,m 

At 

(23a) 

x k+k - 
n,m 

1 (x k 

2 n,m 

+ x k+1 
n,m 

9x k+1 

+ <jf> 4t> 

n,m 

(23b) 


• 3 V 

where x xs the unknown vector to be solved, and (-^) is 

O t 

given by Eqn (14) and (15) . To achieve the formal second- 

3x k 9x k+1 

order accuracy in both space and time, (^ ) and (^£-) are 

computed using one-side difference quotients to replace the 
space derivatives, and alternations between the backward and 
the forward formulas are to be made in Eqn (23a) and (23b) 
for spatial derivatives. On the boundaries of the computa- 
tional region, i.e., body and shock or outer surface due to 
the lack of mesh points outside of the region , alternations 
of the difference quotient can not be applied. However, 
physical boundary conditions are available for these boundaries 
and are taken into account in a two-step fashion consistent 
to the calculations inside the region. Note, the one-side 
difference formulas are only applied to the space derivatives 
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of convection terms with respect to x. However, for the 
stress tensor and the heat flux and mass diffusion vectors, 
centered-dif ference formulas are used. 

It has been recognized that all numerical schemes intro- 
duce a certain amount of numerical dissipation and dispersion 

to the equations being solved and care must be taken to 

21 

distinguish between the numerical and physical dampings 
For example, a first-order scheme could be used for wake 
flow studies if the flow Reynolds number is so low as to 
dominate the flow characteristics. On the other hand, for 
a high Reynolds number blunt body flow, a second-order scheme 
may not be sufficient to suppress the contribution from the 
numerical dissipation. Theoretically, one can always apply 
a fourth or higher-order scheme to the blunt body flow 
problem but the higher-order scheme involves tedious 
arithmetic computations and has not been developed to a 
stage for practical use. Therefore, the only immediate 
remedy for the user is to use very small mesh spacing. 

4 . 2 Numerical Relaxations 

Because of the complexities of the governing equations (14) 
and (15) , it is advantageous to view the time-marching method 
not as a means to solve the mathematical initial-value problem, 
but rather as a relaxation method similar to the one which 
solves the elliptic equations of a boundary-value problem. 

The time-marching method yields the unique steady solution 
no matter what the initial conditions are as long as these 
conditions are compatible and reasonable to the physical 

22 

problem. This point of view was first suggested by Crocco 
and substantiated by later work on both inviscid and viscous 
problems. Interpretation of the time-marching method in 
this way is most suitable to the problem considered in this 
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study, since the transient solution of the problem is of 

no concern. Therefore, the Courant-Friedrich-Lewy (CFL) 

stability condition for a hyperbolic system of equations, 

which provides the largest possible time step size for all 

meshpoints in the computational region for the transient 

problems, will not be followed strictly. Rather it is 

proposed to march the solution in time on each mesh point 

according to its largest possible CFL step size, or to 

enforce the CFL condition by a local value and not by the 
(23) 

global value. This procedure is similar to the variable- 

time method of characteristics whereas the conventional 

procedure corresponds to the fixed-time method of character- 
(24) 

istics. If only the time-asymptotic solution is of 

interest, then the new procedure would reach its limit by a 
smaller number of time steps and the cost of computer time 
would be less than that of the conventional procedure. 

There are additional freedoms available to the users in the 
consideration of time step sizes. If the CFL condition gives 
the relaxation time increment on the basis of the physical 
argument that the propagation of numerical signals should be 
larger than the flow velocity and the speed of sound across 
a mesh spacing, then there exists another relaxation time which 
is valid for chemically reacting flows. The reasoning is this: 
the time step size can not be larger than the chemical 
relaxation time across a mesh spacing in order to maintain 
the stability in the integration of the chemical kinetic 
equations. In the conventional procedure the smaller of 
the two time step sizes is to be used for the whole region. 

But in the new variable-time procedure both the CFL condi- 
tion and chemical relaxation condition are applied locally 
to each mesh point and to different equations. Numerical 
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experiments using this relaxation procedure have been carried 
out for selective cases to be discussed in Section 4.5. 

4 . 3 Sharp-Shock Formulation 

The blunt-body flow field is characterized by the presence 
of gradients of flow properties upstream of the body. At 
high speed and low altitude the gradients are confined in 
a very thin-layer, and can treated mathematically as a 
discontinuity across which the Rankine-Hugoniot relations 
can be applied. Since there is no need to know the struc- 
ture of the shock wave in the computation of the flow field 
downstream of the shock, the shock itself is treated as 
a boundary of the region of interest. The other boundaries 
consist of the axis of symmetry, the body, and a line 
located in the supersonic region connecting the shock and 
the body. The shock wave has been observed experimentally 
as a thin surface so it is assumed that the flow remains 
chemically frozen as it traverses the shock and the 
diffusion in the shock is negligible. In general, a 
correlation can be made between the flow Reynolds number 

5 

and the existence of a sharp shock; namely for Re/Rn >_ 10 
-10 4 the assumption of a sharp shock is very good. Because 
of the high Reynolds number, the flow downstream of the 
shock is largely inviscid with the viscous effects located 
in a thin boundary layer adjacent to the body surface. 

The numerical procedure which relates the shock boundary 
conditions to the finite-difference solution are identical 

(19) 

to the one developed for inviscid flow computations. 

The initial flow field conditions are based on stationary- 
shock results at one side of the computational region and 
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the assumed surface properties on the other. Flow properties 
between these two boundaries are obtained from a linear 
interpolation procedure. The shock location as well as the 
shock speed changes in the subsequent time steps. The shock 
speed, after the first one hundred steps computation, can be 
used to indicate the steadiness of the solution. When the 
ratio of the shock and the free stream speeds becomes less 
than one per cent, the solution can be regarded as the 
steady results. 

To determine the shock speed on various points of the shock 
boundary, a locally intrinsic coordinate system is used, 
while the normal component of the shock speed is used in the 
Rankine-Hugoniot relations. Figure 3 gives the relations 
between the shock intrinsic coordinates and the three possible 
orthogonal coordinates. Details of the matching procedure 
are not different from the discussion in the aforementioned 
references. 

4 . 4 Thick-Shock Formulation 

It has been observed experimentally that the width of the 
gradients of flow properties upstream of a blunt body 

(25) 

increases as the free stream Reynolds number decreases. 

The mechanisms that cause the broadening of a shock wave are 
primarily due to the fact that the ambient air density is 
low and the physical dissipation effects become dominant 
in the flow field. The characteristics of the blunt body 
flow changes drastically as the altitude increases. The 
boundary layer begins to thicken while the shock may still 
be thin, then at higher altitudes the shock width increases 
and merges with the boundary layer. As the altitude further 
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increases there is no distinction between the shock and the 
boundary layer. Thus, the sharp-shock formulation fails to 
provide adequate analysis for flow field predictions. 

The low-Reynolds number blunt body flow problem has been 

studied by many investigators. The continuum approach 

describing the flow field on the basis of the Navier-Stokes 

equations has provided surprisingly accurate results for 

( 26 ) 

Reynolds numbers as low as Re^/Rn =1. ' For slightly 

°° 2 3 

higher Reynolds numbers, say Re^/Rn =10 -10 , the shock 
and boundary layer are not completely merged and the inte- 
gration of equations through the shock as used in Ref. 26 

/ g \ 

and 2 7 can be replaced by the so-called "two- layer*" model. ' 

The computation of the flow field is made downstream from 

the inner-edge of the shock, and the dissipative Rankine- 

Hugoniot relations are used on that edge. Solution can be 

thus obtained in a manner similar to the sharp-shock formu- 
. (14) 

lation. The drawbacks of this method of solution are: 

1) it gives accurate information upstream of the body only 

for certain conditions and the inner-edge of the shock wave 

( 2 7) 

is not physically locatable; 2) it does not clearly 

define the species concentrations on the inner-edge of the 
shock. Ambiguities are introduced to the calculation of species 
on the shock boundary as demonstrated by various investigators. ^ 

The thick-shock formulation adopted in the present study is 
similar to the one for ideal or equilibrium-air-flow calcu- 
lations shown in Ref. 26. The initial conditions are given 
within the computational region which has an outer-boundary 
located far upstream of the body where it is free of any 
disturbance from the body. The free stream properties are 
maintained on this outer-boundary independently of the time 
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steps. The flow properties inside the region are also free 
stream quantities initially and the body is instantaneously 
inserted into the free stream. Note t the use of the outer- 
boundary serves the same purpose as the shock boundary in 
the sharp-shock formulation. The finite-difference method 
is then used to advance the solution in time for mesh points 
inside the computational region, except for the outer-surface. 
The steadiness of the solution is indicated by the negligible 
changes of flow properties between two time steps. 

4 . 5 Verification of the Calculation Procedure 

The basic formulations and the method of solution discussed 
in previous sections have been coded in Fortran V for the 
UNIVAC 1108 system at NASA-MSC. Four versions of the 
program were developed in this study. They are the non- 
reacting viscous thin and thick shock codes, and the 
reacting viscous thin and thick shock codes. These versions 
can be made as options from a general viscous reacting blunt 
body program with some additional programming effort. 
Verification of the basic formulation has been concentrated 
on the non-reacting flows, because more reliable results are 
available for comparison purposes. 

The concept of time-marching according to the local CFL 
step size was first demonstrated in a high Reynolds number 
flow calculation. The sharp-shock formulation was used 
because the steadiness of the solution can be easily judged 
from the magnitude of the shock speed. The free stream 
conditions for this case are shown in Figure 4. A spherical- 
cone is chosen as the body and the downstream outflow boundary 
of the computational region is defined by 9 = 80°. Flow 
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speeds across the outflow line are mostly supersonic. A large 
angle is used here to emphasize the difference between 
solutions which use different procedures of selecting time 
increment. Figure 4 shows the rate of convergence for two 
solutions; one proceeds according to the global CFL time 
increment, whereas the other proceeds according to the local 
CFL time increment. The range of shock speed is given by the 
difference between two curves. It is zero initially and 
should go to zero after many computational steps. With the 
local CFL increment the solution is very close to the steady 
solution after 250 steps. But the other solution indicates 
that more steps are needed before the steady solution is 
reached. Since the global CFL increment is determined most 
likely from the mesh points on the axis of symmetry, these time 
increments are several factors smaller than the CFL increment 
determined from mesh points on the downstream outflow boundary. 
The large difference in the magnitude of the CFL increments 
affects the numerical relaxation to the same degree. Another 
calculation not shown in this report indicates that with a 
smaller computational region, 6 = 60°, the improvement using 
the local CFL increment over the global CFL increment is not 
as great as shown in Figure 3. It is also found that the 
resulting steady solutions are in very close agreement. More 
investigations on the relaxation time increment will be pub- 
lished in a later report. 

The second test case demonstrates the necessity of performing 
flow field computations on the second transformed computational 
plane for a high Reynolds number flow. The body is a 
hyperboloid of 10° asymptotic angle in the free stream 
conditions shown in Figure 5(a) and (b) . Both skin friction 
and heat transfer coefficients were calculated by means of 
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the numerical techniques discussed in Section 2.2 and were 
summarized in Ref. 28. Figure 5 indicates that more than 
100 mesh points in the first computational plane are needed 
across the shock layer, i.e., N=100, in order to predict 
the boundary layer accurately. If the second computational 
plane is used, and with the use of 8=3 and N=45, better 
results can be obtained. The saving of computer cost is 
more than a factor of two because the difference in N used. 
However, it should be pointed out that in the second plane 
the mesh spacing between the shock boundary and the nearest 
mesh line is stretched to a higher degree than the squeezing 
of mesh spacing between the body and its nearest mesh line. 

The stretching of mesh spacing at the shock gives rise to 
difficulties in the marching procedure, and convergence 
of the solution may not be achieved if 8 is too large 
and N is too small. 

It is also of interest to compare other flow properties 
obtained from the non-conservative and con vervative- form 
of governing equations. The governing equations used in 
this work and in Ref. 14 represent the two forms of these 
equations. Figure 6 gives the density and temperature profiles 
on the axis of symmetry for the same free stream conditions. 

The temperatures are quite close, but the densities are very 
different. The density profile obtained from the non-conserva- 
tive form of the equation is unreasonable, especially at the 
body. The fundamental difficulty can be traced back to the 
dependent variables used. Since p, u, v and s, where s is 
the specific entrophy were used in Ref. 14, the boundary 
condition of temperature could not be imposed directly 
upon the solution. Except for this discrepancy, other flow 
quantities are about the same. 
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The last verification case was made by solving the flow field 
structure around a sphere in a low Reynolds number flow. 

The free stream conditions are given in Figure 7 and the 
computational region is given in Figure 1. All the com- 
putations are carried out in the first computational plane, 
since sharp gradients of flow properties disappeared in this 
rarefied gas regime. The transient solution is given in 
Figure 7a. It is observed that at least 1000 time steps 
are required to reach the asymptotic steady solution. The 
overshoot shown in temperature profiles close to the body 
decreases with the increase in time steps. The cause of 
the overshoots is not clear although it is most probably 
from a numerical rather than a physical source. The density 

profile on the axis of symmetry agrees very well with 

( 25 ) 

experimental results . Also shown in Figure 7b are the 

(29) 

results of a Monte Carlo simulation technique . They 

are not as good as the present results. More results for 
a cylinder flow can be found in Ref. 26, 
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5.0 NUMERICAL SOLUTIONS OF THE SHUTTLE CONDITIONS 

A trajectory that corresponds to the long-range orbiter and 
the maximum heating load was used for the flow field compu- 
tation discussed herein. Six representative points in the 
trajectory were selected for the purpose of investigating 
the characteristics of the flow field. The locations of 
the six points are given in the velocity-altitude diagram 
shown in Figure 8. Additional free stream conditions are 
listed in Table 5 for reference. The orbiter trajectory 
covers a wide range of flow regimes that can be categorized 
according to the free stream Reynolds number based on the 
nose radious of 2 ft. Points 0 and 1 are the two highest 
Reynolds numbers flow and are in what is usually called the 
boundary-layer regime. The shock can be treated as a thin 
discontinuity and be treated inviscidly . The viscous 
dissipation is confined in the boundary layer close to the 
nose surface. Points 4 and 5 have the smallest Reynolds 
numbers and are in the rarefied gas regime. A thick shock 
structure appears upstream of the nose and the viscous flow 
extends from the nose to the free stream adjacent to the 
shock. The middle portion of the trajectory is represented 
by points 2 and 3. The shock may be treated as a discon- 
tinuity, however, the boundary-layer is sufficiently thick 
so that interactions between the shock and the boundary 
layer exist. In addition, because of the speed, significant 
departure of chemical equilibrium occurs in the shock layer. 
For the other four points either an equilibrium air or an 
ideal gas model would be the appropriate flow chemistry model 
for both low and high altitude flight, because the speed 
is rather low at the low altitude trajectory while the density 
is low at the high speed and high altitude trajectory. 
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5 . 1 Ideal and Equilibrium Flow Field Solutions 

All the trajectory points except for point 5 were considered 
in the flow field computations using both ideal gas and 
equilibrium air models. A modified set of boundary conditions 
including velocity slip will be used in a future analysis 
to deal with the flow field solution of point 5 exclusively. 

The present analysis is unable to predict physically correct 
results; viz. the skin heat transfer coefficients is larger 
than one, etc. Furthermore, due to the difference in shock 
structure, the sharp shock formulation was used for the 
first four points. Whereas the thick-shock formulation 
was used for point 4. An attempt was made to calculate the 
flow field for point 4 using the sharp-shock formulation 
but the solution could not converge to the steady solution 
due to the small value of the Reynolds number. 

Figure 9a and 9b show the temperature profiles on the axis 
of symmetry obtained for the ideal gas model and the equili- 
brium air model, respectively. Note that the wall temperatures 
are different among the four points in order for the highly- 
cooled wall assumption to be valid. The shock stand-off 
distances are also given in Figure 9a and 9b. The equili- 
brium temperature profile for point 4 is not given because 
of the limitation in the real gas subroutine. 

Figure 10 gives the skin friction coefficients obtained from 
the ideal gas and the equilibrium air analyses. The air 
chemistry only affects the value of Cp, at high altitude. 

It is also found in Figure 11, which presents the skin heat 
transfer coefficients for several points, that the depar- 
ture of chemical equilibrium affects solutions of points 
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2 and 3. It should be pointed out that the highest heat flux 
to the wall is at point 3 computed on the basis of the ideal 
gas and equilibrium air gas models. The nondimensional 
formulas used are 


ns 

I7TT 


GO 00 


(24) 


C H = 


‘n 


[(H ) rn“ (H.J _ J 


W T Y M 


(25) 


where H = e+p/p , and e , are in body intrinsic 

lib XI 

coordinates . 

The results presented in Figure 9 , 10 , and 11 were obtained 
using the local time increment to advance the solution until 
it has reached its asympotic state, and using the second 
computational plane with coordinate squeezing toward the 
body. The mesh was constructed by N=45, M=10, and 8=3. 

The number of time steps needed were K=500-800. The 
execution time of the programs on the UNIVAC 1108 system 
was 45 to 75 minutes for each case. A large number of time 
steps was found to be necessary to bring a steady solution 
to the thick-shock flow field analysis. 


The thermodynamic properties of the ideal gas and the 
equilibrium air are obtained from the following relations, 
p = ( Y _i) pe , T = ( y-1) e/R where y and R are the 
ratio of specific heat and the gas constant, respectively. 
y=1.4 is used in the ideal gas flow field calculations and 
y=y eff is used in the equilibrium air calculation. Y eff 

is determined as a function of internal energy and density 

(30) 

by the following equation : 
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Y eff = R 0 [°- 161 + °* 255 F 1 G 1 + °* 28 ° (1_F 1 ) G 2 + 0,137 F 2 G 3 ] + 1 
where 

a = 0.048 T 1 log 1Q E + 0.032 (1-F ] _) (1-F 2 ) log 1Q E + 0.045 F 2 

and 

E = local internal energy (ergs/gm) X 10 

r = d/d where p = local density, p = sea level density 
o o u 

Also 

E x = 8.50 + 0.357 log 1Q R Q , E 2 = 45.0 R q 0.0157 
and 


AE X = 0.975 


°.° 5 

R o 


AE 0 = 4.0 R 


0.085 


Also i _ 1 

F 1 = [exp (E-E 1 )/AE 1 + l] ' F 2 = [ ex P ( " E + E 2 ) / AE 2 + 1 ] 

G = exp (-E/4.46), G 2 = exp (-E/6.63), and G 3 = exp (-E/25.5) 


5.2 Finite-rate Reacting Flow Field Solutions 

Flow field solutions of trajectory points 2 and 3 are dis- 
cussed in more detail in this section. The effect of the 
body temperature conditions , the surface catalycity , and 
the transport coefficients of the mixture on the flow field 
solution were examined. Comparisons of flow properties 
were also made between the viscous and the inviscid analyses 
and between the different gas models. To simplify the 
investigation, limiting cases of wall conditions are used 
such as an isothermal or adiabatic wall, a fully— catalytic 
or non— catalytic wall. In addition, transport coefficients 
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such as those calculated by the procedure described in 
Section 2.2 or those simply taken from the Sutherland for- 
mula for viscosity coefficients and constant Prandtl and 

Lewis numbers (Pr=0.71 and Le. =1.5) were used, 

1m 

u/y. = 1.5(T/T oo ) 1 * 5 /(T/T m +0.5) 

k = Cpy/Pr , D Zm = k Le Zm /pc p 

The thermodynamic properties of nonreacting gas models are 
obtained from p = (y-l)pe and T = (y-1) e/R/JV^ where 
y = 1.4 for an ideal gas and y = Y e f f for an equilibrium 

air, Y e ff is calculated by a procedure described in refer- 
ence 25. The chemical rate constants used in the reacting 
gas model are given in reference 19. The influence of rate 
constants on the flow properties is not studied in this work, 
since extensive investigations are available in the inviscid 
flow calculations. 

The temperature profiles along the stagnation line and a line 
normal to the wall at S/R^ = 0.7 are shown in figure 12 
for the trajectory point 2 . The temperature profile 
obtained for an ideal gas model exhibits large oscillations, 
especially along the stagnation line where the temperature 
drops quickly to the specified wall temperature. In the 
case of an equilibrium air model, the temperature behind 
the shock is nearly half of that for an ideal gas, and the 
oscillations are relatively small. The finite-rate chemistry 
considered in the shock layer results in a rapidly decreasing, 
but smooth profile toward the wall. This curve resembles 
the result in a fully-viscous shock layer, except for the 
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sharp gradient near the wall. The dash line is the inviscid 
temperature profile obtained from an improved version of 
program reported in reference 19. The inviscid and viscous 
reacting temperatures are very close to each other along 
the stagnation line, but not so close downstream. The 
difference immediately behind the shock is a result of the 
different shock fitting procedures used, however, the 
present one is more accurate. From the shock locations one 
can obtain the velocity displacement thickness or the velocity 
boundary layer thickness. The boundary layer thickens as 
the flow moves downstream and occupies about 10 percent of 
the shock layer. The viscous shock location is affected 
qualitatively by the finite-rate chemistry as the shock in 
inviscid flows. The temperature boundary layer is about 
twice as thick as the velocity boundary layer, and much 
thicker than that obtained from non-reacting calculations. 
Nevertheless, the temperature gradient near the wall is 
still quite high because of the cooled wall. 

The boundary-layer features are also noted in species con- 
centration shown in figure 13. The effect of wall cata- 
lycity is examined for both adiabatic and isothermal walls. 

It is seen that the influence of wall catalycity is restricted 
to ab out one quarter of the stagnation line. For an adiabatic 
and non-catalytic wall, the species profiles of C q2 and 

C „ exhibit very close resemblance to the inviscid species 
NO 

profiles, which are also shown in figure 13a. If the wall 
condition is changed to catalytic, the influence of wall 
catalycity is extended to one- third of the stagnation line. 
Since the flow properties are affected to different degrees 
by the wall conditions, it is therefore difficult to obtain 
appropriate edge conditions for boundary-layer analyses. 
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Another flowfield calculation was made using multicomponent 
transport coefficients for a high-cooled, catalytic wall. 
Figures 14 and 15, respectively, present the temperature 
profiles across the shock layer at four points on the body 
and species profiles at two points on the body. The flow- 
field results are nearly indistinguishable from each other 
as comparison is made between figures 12 and 14 and between 
13 and 15. In terms of computation time, however, the 
multicomponent results requires four times as much as that 
for results obtained from using simple transport coefficients. 
The computation time per mesh point is about one-tenth of 
a second for viscous calculations and about one-fifteenth 
of a second for inviscid calculations. It is estimated for 
a reacting flow calculation using seven species and 19 
reactions. The total time for a converged solution is 45 
minutes for a mesh of 10x15, or 10 points along the body 
and 15 points across the shock layer, with 250 time steps. 

More computation time can be saved for multicomponent solution 
by matching the inviscid and viscous solutions outside of 
the boundary layer. 

Because of the cost to execute the program, most of the 
results are obtained from a mesh of 10x15. The resultant 
spacing near the wall is obvisouly too coarse to predict 
accurate heating rates, but is sufficiently small to resolve 
the fine detail of the flow property profiles. The overall 
accuracy on the flowfield solution is checked by comparing 
the total specific energy with its free stream value. The 
dissipations of total energy are presented in figure 16 at 
three body points. All three curves decrease rapidly in 
proportion to the distance from the wall. The level of 
numerical dissipation is also given for an inviscid flow 
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calculation of same free stream condition. The accuracy of 
the computed flowfield results is then checked by using a 
finer mesh spacing. Difference in flow properties are 
within 5 percent between the solutions obtained from 10x15 
and 10x30. However, the mesh spacing is critical to the 
calculation of skin friction and heat transfer coefficient. 

It is shown in figure 21 that the C„ obtained using a mesh 

H. 

of 10x15 is almost one order of magnitude less than the 

( 2 ) 

existing results. Improvement can be achieved by re- 

ducing the mesh spacing near the wall. It is found that 
the C„ values begin to level off with a mesh of 10x30 and 
3=2, where the mesh spacing next to the wall is approxi- 
mately one-hundredth of the shock layer thickness. More 
discussion on the calculations with a stretched coordinate 
can be found in reference 14. It should be pointed out that 
the multicomponent C„ is higher than the simple C„ by a 

n n 

factor of 2 or more. This is caused by the higher coefficients 
of mass diffusion and thermal conductivity of the multi- 
component theory. 

The flowfield results of the trajectory point 3 are presented 
in figures 17 through 20. The ideal gas temperature profile 
has larger oscillations and the reacting temperature drops 
more quickly behind the shock. Also, the boundary-layers 
are thicker than those of the trajectory point 2. There 
is roughly a factor of two increase for both velocity and 
temperature boundary-layer thickness. Due to a much higher 
speed, the dissociation inside the shock layer is more 
complete and ionization is considered in the calculation. 

But the boundary-layer features are still observed in the 
flow property profiles. Notice that the peak of C n in 

'-'O 
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figure 17 is located approximately in the middle of the 
shock layer for three combinations of wall conditions. The 
interaction between the boundary layer and the shock is 
markedly stronger than that for the trajectory point 2. 

Figures 19 and 20 give the temperature and species profiles 
across the shock layer at several points on the body. Very 
little difference is found between the results obtained from 
the simple transport theory and the multicomponent theory. 
Finally, the heating rates are computed and presented in 
figure 21. Note, the agreement with existing data is not 
quite good. The calculations of C„ have been made on 
meshes of 10x20, 3 = 1 and 10x35; g = 2 . But the value 

of C„ is still lower than the boundary-layer result. 

Some of the results obtained for the trajectory point 4 are 
shown in figures 22 and 23. The assmuption of a thin shock 
was used in the reacting air calculation with simple trans- 
port properties. The shock layer is fully viscous, as the 
wall effect extends to the shock and the temperature profile 
no longer exhibits markedly decreasing magnitude down- 
stream of the shock. In contrast to the temperature pro- 
file obtained for the trajectory points 2 and 3, the tempera- 
ture gradients at the wall are also smaller. Another 
interesting result is that the shock wave is located a little 
farther than that of the trajectory point 3. The free stream 
Mach number are, respectively 27- and 29.5 for points 3 and 
4. The validity of applying the thin-shock program to the 
trajectory point, which has a Reynolds number of 900 based 
on a nose radius of 2 ft, was checked by using the thick-shock 
program. The dash line presented in figure 22 are the ideal 
gas results for the trajectory point 4. Significant 
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differences are observed between these two sets of calcula- 
tions. In addition to the effects of different gas models 
used, the thick-shock calculation was made for Re = 450 
or for the same free stream condition but one foot nose 
radius. The computation encountered difficulty in reaching 
the steady state for Re^ = 900 or higher Reynolds numbers. 

It seems that a certain amount of physical dissipation is 
needed to damp the oscillations in the flowfiled results for 
a chosen mesh size. The mesh used in the present calcula- 
tion is N = 20 and M = 10 and the computational region is 
defined by S/R^^ = 0.6 N/Rj^ = 1 • The computer time is 
about 85 minutes on the MSC-UNIVAC 1108 for K = 1500 . 

Since the reacting air calculation would require several 
fold more computer time, such a calculation has not been 
attempted. 

Figure 23 shows the species profiles on two stations across 
the shock layer obtained from the thin-shock reacting program. 
The level of dissociation and ionization is approximately the 
same as that of the trajectory point 3. However, the shape 
of the profiles are quite different as revealed by a com- 
parison between Figures 23 and 20, where the wall conditions 
are the same for both cases. 

Finally, in Figure 24 the best rates are summarized for the 
trajectory points considered in this study. 

The results for points 0 and 5 are not presented because 
the speed is either too low or the altitudes too high for 
these conditions such that the departure of chemical equili- 
brium is negligible in the shock layer. 
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The value of C„ tends to agree with the results obtained 
H , 

from the boundary-layer or the thin shock layer theories, 
although our value is consistently lower for all the points. 
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6.0 CONCLUSION 

A time-marching finite-difference method has been developed 
to solve for the flow field in the nose region of an orbiter. 
During high altitude flight the thick-shock formulation predicts 
results that match with the free-molecular limit, whereas at 
lower altitude the thin-shock formulation with the aid of 
coordinate transformations provides flow field solutions that 
agree quite closely with the usual boundary- layer theory. 

The present method is particularly useful to analyze the 
flow field for the medium altitude when appreciable chemical 
nonequilibrium exists. Since an "exact" form of the governing 
equations is utilized and only numerical approximations are 
made in the method of solution, the solutions obtained are 
considered to be more general and satisfactory than most 
results obtained from nonequilibrium thin shock-layer theory. 

The shortcomings of existing theories; i.e., the need of 
accurate edge conditions for reacting boundary- layer analysis, 
and the use of a "two-layers" model and of simple transport 
coefficients analysis are absent from the present solution. 

Two major assumptions are employed in the present, as well 
as, in the previous analysis; namely that the vibrational 
equilibrium is maintained and mass diffusion depends only 
on the density gradient. These assumptions are made in 
order for the computers presently available to manage the 
computations within a reasonable time. The cost of the 
flow field computation is noted to increase considerably 
with the use of a more elaborate theory of the chemical 
kinetics and transport properties. The application of the 
time-marching method to more practical problems rests 
entirely on the further improvement of its efficiency. It 
is felt that more studies should be directed toward the 
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numerical experimentation that will provide new techniques 
with the potential to expedite the convergence to the time- 
asympotic solution. Also, the accuracy of the results thus 
obtained should be examined against the experimental data 
which are not currently available. 
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APPENDIX I 

TEMPERATURE COEFFICIENTS FOR THERMODYNAMIC FUNCTIONS 

The first line indicates the range of temperatures for 
which the following polynomials are fitted. The first two 
numbers give the first temperature interval, the last two 
give the second temperature interval. The second line 
refers to the name of species that is followed by the mole- 
cular weight of the species, M^the mass concentration in the 
free stream, C. , and the molar enthalpy of the species at 

oo 

temperature °K, H.. The third to the fifth lines are two 

O 

sets of temperature coefficients for the temperature inter- 
vals stated in the first line. Each set of the temperature 
coefficients contains seven numbers, A^ (i=l,...7). 

The input to the "Fortran IV program for calculation of 
thermodynamic data" ^ 16) is also given for reference. 
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TEMPERATURE COEFFICIENTS ffOR THERMODYNAMIC FUNCTIONS 
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APPENDIX II 

CROSS SECTIONS FOR CALCULATION OF TRANSPORT PROPERTIES 

The cross section data are obtained from Ref. 18 for the 
neutral species of air. Each pair of the interactions is 
listed in the first line, also shown is the number of 
temperatures followed. For each temperature, ' 

A^ m and are given. The subroutines responsible for 

the calculation of transport properties of the air mixture 
are written on the basis of these developed in Ref. 18. 
Several important modifications and numerous changes in 
coding have been made to the original version in order to 
suit the flowfield program. 
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APPENDIX III 

CHEMICAL REACTIONS AND RATE CONSTANT 

The calculation of the net rate of chemical production, 
for 1-species is made in one subroutine which was 

written previously for the inviscid flowfield calculations. 
The rate constants given are taken from Ref. 31. 
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TABLE 2. NUMERICAL METHODS FOR ANALYZING VISCOUS REACTING FLOW FIELD 
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ORTHOGONAL COORDINATE SYSTEM AND ITS SPECIAL FORMS 
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METRIC COEFFICIENTS AND THEIR DERIVATIVES IN THE COMPUTATIONAL PLANE 







TABLE 5. FREE STREAM CONDITIONS FOR THE ORBITER TRAJECTORY POINTS INVESTIGATED 

IN THIS STUDY 
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TABLE 6 . DESCRIPTIONS OF CASES CALCULATED FOR AN ORBITER 
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